Propofol Dosage Analysis IU Health Ball Memorial Hospital
Jasmin Cooper
2025-10-20
#load datasets and libraries
library(rsample)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.4
## ✔ ggplot2 3.4.4 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(naniar)
library(corrplot)## corrplot 0.92 loadedPropdata <- read_csv("propofoldata.csv")## Rows: 72 Columns: 44
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (23): Org_Name, Height, IBW, Triglyc_Max, CrCl, HR_Range_low, Min_Rass, ...
## dbl (21): Age, Sex, Weight, COVID, IVDA_Hx, Home_Opiod, Home_Benzo, Opiate_T...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(Propdata)## # A tibble: 6 × 44
## Org_Name Age Sex Height Weight IBW COVID IVDA_Hx Home_Opiod Home_Benzo
## <chr> <dbl> <dbl> <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 IU Health… 43 0 154.9 118. 47.7… 0 0 1 0
## 2 IU Health… 65 0 149.9 73.9 43.2… 0 0 0 0
## 3 IU Health… 64 0 180.3 128. 70.7… 1 0 0 0
## 4 IU Health… 67 0 157.5 155. 50.1… 0 0 0 0
## 5 IU Health… 63 0 167.6 79.4 59.2… 1 0 0 0
## 6 IU Health… 37 0 172.7 164 63.8… 0 0 1 1
## # ℹ 34 more variables: Opiate_Tol <dbl>, Triglyc_Max <chr>, Creatinine <dbl>,
## # CrCl <chr>, Map_Range <dbl>, HR_Range_low <chr>, Total_days_prop <dbl>,
## # Dur_prop_minutes <dbl>, Per_prop_over50 <dbl>, Max_Dose_Prop <dbl>,
## # Min_Rass <chr>, Max_Rass <chr>, Dex_ord <chr>, Prop_dose_when_Dex <chr>,
## # Dex_rate_when_prop <chr>, Prop_dec_afterDex <chr>, Fent_before_prop <chr>,
## # Prop_dose_when_fent <chr>, Fent_max_when_prop <chr>,
## # Prop_dec_after_fent <chr>, Numb_fent_bolus <chr>, Mode_CPOT <chr>, …#show(Propdata)
#print(Propdata$Opiate_Tol)Creating Analysis Dataset
IBW, COVID, IVDA_Hx, Triglyc_Max, Creatinine, CrCl, Map_Range, and HR_Range_low will all be dropped from the analysis dataset since they are not relevant to the analysis being completed.
#creating analysis dataset
df_eda <- select(Propdata, -c(IBW, COVID, IVDA_Hx,Triglyc_Max,Creatinine,CrCl,Map_Range,HR_Range_low))
df_eda[df_eda=="x"] <- NA
df_eda[df_eda=="X"] <- NAx/X that represent NA/missing values is replaced with NA to be able to visualize and replace NA values.
#visualizing and appropriately replacing missing variables
vis_miss(df_eda[,1:12])
df_eda$Height <- ifelse(is.na(df_eda$Height), "Missing", df_eda$Height)
vis_miss(df_eda[,13:25])
df_eda$Min_Rass <- ifelse(is.na(df_eda$Min_Rass), "NotApp", df_eda$Min_Rass)
df_eda$Max_Rass <- ifelse(is.na(df_eda$Max_Rass), "NotApp", df_eda$Max_Rass)
df_eda$Dex_ord <- ifelse(is.na(df_eda$Dex_ord), "NotApp", df_eda$Dex_ord)
df_eda$Prop_dose_when_Dex <- ifelse(is.na(df_eda$Prop_dose_when_Dex), "NotApp", df_eda$Prop_dose_when_Dex)
df_eda$Dex_rate_when_prop <- ifelse(is.na(df_eda$Dex_rate_when_prop), "NotApp", df_eda$Dex_rate_when_prop)
df_eda$Prop_dec_afterDex <- ifelse(is.na(df_eda$Prop_dec_afterDex), "NotApp", df_eda$Prop_dec_afterDex)
df_eda$Fent_before_prop <- ifelse(is.na(df_eda$Fent_before_prop), "NotApp", df_eda$Fent_before_prop)
df_eda$Prop_dose_when_fent <- ifelse(is.na(df_eda$Prop_dose_when_fent), "NotApp", df_eda$Prop_dose_when_fent)
df_eda$Fent_max_when_prop <- ifelse(is.na(df_eda$Fent_max_when_prop), "NotApp", df_eda$Fent_max_when_prop)
df_eda$Prop_dec_after_fent <- ifelse(is.na(df_eda$Prop_dec_after_fent), "NotApp", df_eda$Prop_dec_after_fent)
df_eda$Numb_fent_bolus <- ifelse(is.na(df_eda$Numb_fent_bolus), "NotApp", df_eda$Numb_fent_bolus)
df_eda$Mode_CPOT <- ifelse(is.na(df_eda$Mode_CPOT), "NotApp", df_eda$Mode_CPOT)
df_eda$Was_ket_ord <- ifelse(is.na(df_eda$Was_ket_ord), "NotApp", df_eda$Was_ket_ord)
vis_miss(df_eda[,26:36])
df_eda$Prop_dose_when_ket <- ifelse(is.na(df_eda$Prop_dose_when_ket), "NotApp", df_eda$Prop_dose_when_ket)
df_eda$Ket_max_when_prop <- ifelse(is.na(df_eda$Ket_max_when_prop), "NotApp", df_eda$Ket_max_when_prop)
df_eda$Was_prop_dec <- ifelse(is.na(df_eda$Was_prop_dec), "NotApp", df_eda$Was_prop_dec)
df_eda$Con_Prop_NMBA <- ifelse(is.na(df_eda$Con_Prop_NMBA), "NotApp", df_eda$Con_Prop_NMBA)
vis_miss(df_eda)
View(df_eda)Exploratory Data Analysis
Visualizing categorical variables
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=Org_Name))
Most patients are from the IU Health Ball Memorial Hospital. This is helpful since this is where we are concentrating our analysis.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Sex), fill=factor(Sex)))
There are more males than females.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Home_Opiod), fill=factor(Home_Opiod)))
Most patients do not have home opiod usage.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Home_Benzo), fill=factor(Home_Benzo)))
Most patients do not have home benzo usage.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Opiate_Tol), fill=factor(Opiate_Tol)))
Most patients have no tolerance for Opiate.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Min_Rass), fill=factor(Min_Rass)))
Min RASS documented when propofol rate >50 mcg/kg/min. “The Richmond Agitation Sedation Scale (RASS) is an instrument designed to assess the level of alertness and agitated behavior in critically-ill patients.” (1) Most patients minimum RASS when propofol rate >50 mcg/kg/min was that they were unarousable. Link to learn more about the scale: https://www.physio-pedia.com/Richmond_Agitation-Sedation_Scale_(RASS)
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Max_Rass), fill=factor(Max_Rass)))
Max RASS documented when propofol rate >50 mcg/kg/min. “The Richmond Agitation Sedation Scale (RASS) is an instrument designed to assess the level of alertness and agitated behavior in critically-ill patients.” (1) Most patients maximum RASS when propofol rate >50 mcg/kg/min was that they would pull at or attempt to remove tubes alongside generally aggressive behavior. Link to learn more about the scale: https://www.physio-pedia.com/Richmond_Agitation-Sedation_Scale_(RASS)
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Dex_ord), fill=factor(Dex_ord)))
Dexmedetomidine was ordered before propofol rate reached >50 mcg/kg/min at any point. Most patients were not ordered dexmedetomidine.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Prop_dec_afterDex), fill=factor(Prop_dec_afterDex)))
Propofol was decreased to <50 mcg/kg/min after dexmedetomidine was added (did not go back above 50). Since most patients were not ordered dexmedetomidine, this is not applicable to most patients.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Fent_before_prop), fill=factor(Fent_before_prop)))
Fentanyl was ordered before propofol rate was >50 mcg/kg/min. Most patients did have fentanyl ordered for them before the propofol rate was >50 mcg/kg/min.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Prop_dec_after_fent), fill=factor(Prop_dec_after_fent)))
Propofol was decreased to <50 mcg/kg/min after fentanyl was added. This was not applicable to most patients.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Mode_CPOT), fill=factor(Mode_CPOT)))
Mode CPOT when propofol rate >50 mcg/kg/min . “The CPOT is a behavioural assessment pain scale for patients unable to verbalise pain.” (2) Most patients had a mode CPOT of 0, which means that overall, the patient was relaxed, had no muscle tension, and had an absence of movement. Link to learn more about the scale: https://ccs-sth.org/resources/Documents/Sedation%20Analgesia%20Delirium%20in%20CC/CCS%20STH%20Guideline%20template%20CPOT.pdf
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Was_ket_ord), fill=factor(Was_ket_ord)))
Most patients did not have ketamine prescribed to them before propofol rate was >50 mcg/kg/min and out of the ones who did, most did not have it ordered.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Was_prop_dec), fill=factor(Was_prop_dec)))
Since most patients did not have ketamine prescribed to them, it was not applicable to most patients if propofol decreased to <50 mcg/kg/min after ketamine was added or the dose was increased.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Rec_pain_med), fill=factor(Rec_pain_med)))
Most patients did not receive any scheduled pain medications when propofol >50 mcg/kg/min.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(On_benzo), fill=factor(On_benzo)))
Most patients were not on scheduled benzo.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Con_Prop_Vas), fill=factor(Con_Prop_Vas)))
More patients did not receive concurrent propofol and vasopressor when propofol >50 mcg/kg/min than those who did.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Con_Prop_Dex), fill=factor(Con_Prop_Dex)))
More patients did not receive concurrent propofol and dexmedetomidine than those who did.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Con_Prop_Fen), fill=factor(Con_Prop_Fen)))
More patients did receive concurrent propofol and fentanyl infusion than those who did not.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Con_Prop_Ket), fill=factor(Con_Prop_Ket)))
Most patients did not receive concurrent propofol and ketamine infusion.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Con_Prop_Ben), fill=factor(Con_Prop_Ben)))
Most patients did not receive concurrent propofol and benzodiazepine infusion.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=factor(Con_Prop_NMBA), fill=factor(Con_Prop_NMBA)))
More patients did not receive concurrent propofol and NMBA infusion than those who did.
Visualizing covariation between two categorical variables
Creating initialisms of Org_Name values to make bar graphs easier to interpret.
df_eda$Org_Name[df_eda$Org_Name=="IU Health Arnett Hospital"] <- 'IU HAH'df_eda$Org_Name[df_eda$Org_Name=="IU Health Ball Memorial Hospital"] <- 'IU HBMH'df_eda$Org_Name[df_eda$Org_Name=="IU Health Bloomington Hospital"] <- 'IU HBloomH'df_eda$Org_Name[df_eda$Org_Name=="IU Health Methodist Hospital"] <- 'IU HMH'df_eda$Org_Name[df_eda$Org_Name=="IU Health University Hospital"] <- 'IU HUH'Initialisms for hospital names:
- IU Health Arnett Hospital <- ‘IU HAH’
- IU Health Ball Memorial Hospital <- ‘IU HBMH’
- IU Health Bloomington Hospital <- ‘IU HBloomH’
- IU Health Methodist Hospital <- ‘IU HMH’
- IU Health University Hospital <- ‘IU HUH’
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Home_Opiod)), position = "dodge")
The organization with the most patients that have home opiod usage is IU HBMH.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Home_Benzo)), position = "dodge")
The organization that has the most patients with home benzo usage is IU HBloomH.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Opiate_Tol)), position = "dodge")
IU HBMH has the most patients with high tolerance of opiates.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Prop_dec_afterDex)), position = "dodge")
Only IU HBMH and IU HBloomH had propofol decreases to <50 mcg/kg/min after dexmedetomidine was added.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Fent_before_prop)), position = "dodge")
IU HBMH had the highest number of fentanyl orders before propofol rate was >50 mcg/kg/min.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Prop_dec_after_fent)), position = "dodge")
Only IU HBloomH and IU HMH had propofol decreases to <50 mcg/kg/min after fentanyl was added.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Was_ket_ord)), position = "dodge")
IU HMH had the highest number of ketamine orders before propofol rate was >50 mcg/kg/min.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Was_prop_dec)), position = "dodge")
Propofol was only decreased to <50 mcg/kg/min after ketamine was added or dose increased at IU HBMH.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Rec_pain_med)), position = "dodge")
Only patients at IU HMH and IU HUH received scheduled pain medications when propofol >50.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Vas)), position = "dodge")
Patients received concurrent propofol and vasopressor when propofol >50 mcg/kg/min at IU HBMH, IU HMH, and IU HBloomH.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Dex)), position = "dodge")
Only IU HUH did not have patients who received concurrent propofol and dexmedetomidine.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Fen)), position = "dodge")
Concurrent propofol and fentanyl infusion is a popular choice at all organizations except IU HUH.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Ket)), position = "dodge")
On the contrary, concurrent propofol and ketamine infusion was not a popular choice at the organizations.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_Ben)), position = "dodge")
IU HMH had the most patients that received concurrent propofol and benzodiazepine infusion.
ggplot(data=df_eda)+
geom_bar(mapping = aes(x=Org_Name, fill=factor(Con_Prop_NMBA)), position = "dodge")
IU HBMH had the most patients that received concurrent propofol and NMBA infusion.
Visualizing continuous variables
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = Age), binwidth = 3)
Most patients are between 50 and 70 years old.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = as.numeric(Height)), binwidth = 10)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).
Most patients are around 180cm in height (roughly 5 ft 11in).
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = Weight), binwidth = 5)
Weights are pretty evenly distributed with a few obvious outliers.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = Total_days_prop), binwidth = 1)
Most patients were on propofol for 15 days or less, with a few outliers towards 20 days and 30 days.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = Dur_prop_minutes), binwidth = 2000)
Most patients were on propofol for less than 10000 minutes (which is roughly less than a week) with a few higher outliers.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = Per_prop_over50), binwidth = 5)
The % of time of propofol above 50 mcg/kg/min is relatively spread across patients.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = Max_Dose_Prop), binwidth = 3)
The max dose of propofol in mcg/kg/min greatly varies across patients.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = as.numeric(Prop_dose_when_Dex)), binwidth = 3)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion## Warning: Removed 52 rows containing non-finite values (`stat_bin()`).
This was not applicable to most patients, but for the patients where it was applicable, the propofol dose in mcg/kg/min when dexmedetomidine was added varied across patients.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = as.numeric(Dex_rate_when_prop)), binwidth = 0.25)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion## Warning: Removed 52 rows containing non-finite values (`stat_bin()`).
This was not applicable for most patients, but for the patients where it was applicable, the dexmedetomidine max rate in mcg/kg/hr when propofol >50 mcg/kg/min was usually around 0.5 mcg/kg/hr.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = as.numeric(Prop_dose_when_fent)), binwidth = 0.25)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion## Warning: Removed 59 rows containing non-finite values (`stat_bin()`).
This was not applicable to most patients, but for the patients where it was applicable, the propofol dose in mcg/kg/min when fentanyl infusion was added was varied.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = as.numeric(Fent_max_when_prop)), binwidth = 5)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion## Warning: Removed 22 rows containing non-finite values (`stat_bin()`).
The fentanyl max rate in mcg/hr when propofol >50 mcg/kg/min was between 100-150 mcg/hr for most patients.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = as.numeric(Numb_fent_bolus)), binwidth = 5)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion## Warning: Removed 13 rows containing non-finite values (`stat_bin()`).
Most patients did not receive fentanyl bolus when propofol rate >50 mcg/kg/min.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = as.numeric(Prop_dose_when_ket)), binwidth = 3)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion## Warning: Removed 65 rows containing non-finite values (`stat_bin()`).
This was not applicable for most patients, but for the patients where it was applicable, the propofol dose in mcg/kg/min when fentanyl infusion was added was usually 50 mcg/kg/min or above.
ggplot(data=df_eda)+
geom_histogram(mapping = aes(x = as.numeric(Ket_max_when_prop)), binwidth = 0.25)## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning in FUN(X[[i]], ...): NAs introduced by coercion## Warning: Removed 61 rows containing non-finite values (`stat_bin()`).
This was not applicable for most patients, but for the patients where it was applicable, the ketamine max rate in mcg/kg/hr when propofol >50 mcg/kg/min was spread between 0.0 to 2.0 mcg/kg/hr.
Visualizing covariation between a categorical variable and a continuous variable
The quantitative response variable I have chosen is Per_prop_over50, which is the % of time of propofol above 50 mcg/kg/min. Based off evaluating the EDA results, the selected covariates (explanatory variables) are: Home_Opiod, Home_Benzo, Opiate_Tol, Prop_dec_afterDex, Fent_before_prop, Prop_dec_after_fent, Was_prop_dec, Rec_pain_med, Con_Prop_Vas, Con_Prop_Dex, Con_Prop_Fen, Con_Prop_Ket, Con_Prop_Ben, and Con_Prop_NMBA. Below are some additional visualizations showing covariation between per_prop_over50 and the selected covariates.
ggplot(data = df_eda, mapping = aes(x = factor(Org_Name), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Home_Opiod), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Home_Benzo), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Opiate_Tol), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Min_Rass), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Max_Rass), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Dex_ord), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Prop_dec_afterDex), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Fent_before_prop), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Mode_CPOT), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Was_ket_ord), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Was_prop_dec), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Rec_pain_med), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(On_benzo), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Vas), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Dex), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Fen), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Ket), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_Ben), y = Per_prop_over50)) + geom_boxplot()
ggplot(data = df_eda, mapping = aes(x = factor(Con_Prop_NMBA), y = Per_prop_over50)) + geom_boxplot()
Fitting GLM by selecting a quantitative response variable column and a set of covariates
The quantitative response variable I have chosen is Per_prop_over50, which is the % of time of propofol above 50 mcg/kg/min. The selected covariates (explanatory variables) are: Home_Opiod, Home_Benzo, Opiate_Tol, Prop_dec_afterDex, Fent_before_prop, Prop_dec_after_fent, Was_prop_dec, Rec_pain_med, Con_Prop_Vas, Con_Prop_Dex, Con_Prop_Fen, Con_Prop_Ket, Con_Prop_Ben, and Con_Prop_NMBA.
glm_mod1 <- glm(Per_prop_over50~ Home_Opiod+Opiate_Tol+Dex_ord+Prop_dose_when_Dex+Fent_before_prop+Prop_dec_after_fent+Mode_CPOT+Was_ket_ord+Rec_pain_med+On_benzo+Con_Prop_Vas+Con_Prop_Dex+Con_Prop_Ket+Con_Prop_Ben+Con_Prop_NMBA ,data=df_eda,family="gaussian")
summary(glm_mod1)##
## Call:
## glm(formula = Per_prop_over50 ~ Home_Opiod + Opiate_Tol + Dex_ord +
## Prop_dose_when_Dex + Fent_before_prop + Prop_dec_after_fent +
## Mode_CPOT + Was_ket_ord + Rec_pain_med + On_benzo + Con_Prop_Vas +
## Con_Prop_Dex + Con_Prop_Ket + Con_Prop_Ben + Con_Prop_NMBA,
## family = "gaussian", data = df_eda)
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.0772 26.0465 0.425 0.67366
## Home_Opiod 11.4687 10.0100 1.146 0.26097
## Opiate_Tol 3.7031 6.8949 0.537 0.59517
## Dex_ord1 17.4662 41.1774 0.424 0.67447
## Dex_ordNotApp 25.4739 40.2750 0.632 0.53185
## Prop_dose_when_Dex10 39.3687 34.3706 1.145 0.26109
## Prop_dose_when_Dex100 29.1390 23.9251 1.218 0.23275
## Prop_dose_when_Dex21.2404 28.3292 28.2648 1.002 0.32423
## Prop_dose_when_Dex40 -11.1913 23.8410 -0.469 0.64217
## Prop_dose_when_Dex50 22.6597 24.5651 0.922 0.36366
## Prop_dose_when_Dex50.0002 7.8072 23.2330 0.336 0.73918
## Prop_dose_when_Dex51.2821 61.9114 57.6052 1.075 0.29105
## Prop_dose_when_Dex52 31.9668 36.3055 0.880 0.38559
## Prop_dose_when_Dex53.8642 45.4851 30.6342 1.485 0.14803
## Prop_dose_when_Dex54.1463 32.0240 36.1578 0.886 0.38284
## Prop_dose_when_Dex59.991 27.6077 24.5027 1.127 0.26879
## Prop_dose_when_Dex60 21.7064 29.6990 0.731 0.47052
## Prop_dose_when_Dex65 9.7494 23.9251 0.407 0.68654
## Prop_dose_when_Dex70 2.7771 21.2857 0.130 0.89707
## Prop_dose_when_Dex80 11.2612 34.7263 0.324 0.74797
## Prop_dose_when_Dex89.8261 26.3817 25.4026 1.039 0.30732
## Prop_dose_when_DexNotApp -14.8992 43.0625 -0.346 0.73177
## Fent_before_prop1 -5.6738 14.8100 -0.383 0.70434
## Fent_before_propNotApp 10.1541 15.8590 0.640 0.52686
## Prop_dec_after_fent1 -19.2199 17.0755 -1.126 0.26927
## Prop_dec_after_fentNotApp 8.4598 16.7658 0.505 0.61753
## Mode_CPOT1 3.4145 24.2573 0.141 0.88900
## Mode_CPOT2 11.4906 12.1827 0.943 0.35312
## Mode_CPOT3 14.5515 11.0369 1.318 0.19733
## Mode_CPOT4 0.8971 12.7054 0.071 0.94418
## Mode_CPOT5 24.4714 9.8122 2.494 0.01837 *
## Mode_CPOT6 37.9780 13.7865 2.755 0.00989 **
## Mode_CPOT7 NA NA NA NA
## Mode_CPOT8 -14.0736 26.9472 -0.522 0.60532
## Mode_CPOTNotApp -3.6458 17.8405 -0.204 0.83945
## Was_ket_ord1 -12.6255 23.6656 -0.533 0.59762
## Was_ket_ordNotApp -23.6351 19.8455 -1.191 0.24301
## Rec_pain_med -14.5930 13.7012 -1.065 0.29533
## On_benzo -8.1254 18.2073 -0.446 0.65861
## Con_Prop_Vas 4.5572 8.1382 0.560 0.57965
## Con_Prop_Dex NA NA NA NA
## Con_Prop_Ket NA NA NA NA
## Con_Prop_Ben -17.5800 9.6445 -1.823 0.07831 .
## Con_Prop_NMBA1 10.8414 6.6535 1.629 0.11368
## Con_Prop_NMBANotApp 8.5955 24.0118 0.358 0.72287
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 287.5668)
##
## Null deviance: 31735 on 71 degrees of freedom
## Residual deviance: 8627 on 30 degrees of freedom
## AIC: 634.92
##
## Number of Fisher Scoring iterations: 2Creating new model with the removal of home_opiod, opiate_tol, dex_ord, fent_before_prop, on_benzo, con_prop_vas, con_prop_dex, and con_prop_ket due to insignificance.
glm_mod2 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Prop_dec_after_fent+Mode_CPOT+Was_ket_ord+Rec_pain_med+Con_Prop_Ben+Con_Prop_NMBA ,data=df_eda,family="gaussian")
summary(glm_mod2)##
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Prop_dec_after_fent +
## Mode_CPOT + Was_ket_ord + Rec_pain_med + Con_Prop_Ben + Con_Prop_NMBA,
## family = "gaussian", data = df_eda)
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.464 19.087 0.181 0.85694
## Prop_dose_when_Dex10 39.541 31.496 1.255 0.21700
## Prop_dose_when_Dex100 30.581 22.727 1.346 0.18641
## Prop_dose_when_Dex21.2404 39.896 23.530 1.696 0.09815 .
## Prop_dose_when_Dex40 10.636 22.727 0.468 0.64247
## Prop_dose_when_Dex50 9.453 22.667 0.417 0.67900
## Prop_dose_when_Dex50.0002 16.457 22.667 0.726 0.47228
## Prop_dose_when_Dex51.2821 63.853 36.876 1.732 0.09146 .
## Prop_dose_when_Dex52 54.704 32.230 1.697 0.09782 .
## Prop_dose_when_Dex53.8642 46.008 25.953 1.773 0.08429 .
## Prop_dose_when_Dex54.1463 59.457 27.126 2.192 0.03459 *
## Prop_dose_when_Dex59.991 34.877 22.157 1.574 0.12376
## Prop_dose_when_Dex60 39.255 24.206 1.622 0.11314
## Prop_dose_when_Dex65 11.191 22.727 0.492 0.62525
## Prop_dose_when_Dex70 18.863 19.250 0.980 0.33334
## Prop_dose_when_Dex80 32.804 30.175 1.087 0.28382
## Prop_dose_when_Dex89.8261 28.054 22.157 1.266 0.21317
## Prop_dose_when_DexNotApp 20.246 14.433 1.403 0.16880
## Prop_dec_after_fent1 -21.775 15.536 -1.402 0.16917
## Prop_dec_after_fentNotApp -1.862 10.634 -0.175 0.86194
## Mode_CPOT1 -4.981 14.812 -0.336 0.73849
## Mode_CPOT2 3.884 11.433 0.340 0.73595
## Mode_CPOT3 9.856 9.335 1.056 0.29774
## Mode_CPOT4 10.898 10.277 1.060 0.29565
## Mode_CPOT5 28.093 8.862 3.170 0.00301 **
## Mode_CPOT6 32.399 13.153 2.463 0.01841 *
## Mode_CPOT7 NA NA NA NA
## Mode_CPOT8 -4.369 23.829 -0.183 0.85549
## Mode_CPOTNotApp -11.875 17.618 -0.674 0.50436
## Was_ket_ord1 5.466 19.175 0.285 0.77714
## Was_ket_ordNotApp -9.113 14.031 -0.650 0.51992
## Rec_pain_med -6.995 12.754 -0.548 0.58661
## Con_Prop_Ben -13.586 8.041 -1.690 0.09931 .
## Con_Prop_NMBA1 8.191 5.391 1.520 0.13691
## Con_Prop_NMBANotApp 6.858 22.667 0.303 0.76389
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 295.4365)
##
## Null deviance: 31735 on 71 degrees of freedom
## Residual deviance: 11227 on 38 degrees of freedom
## AIC: 637.88
##
## Number of Fisher Scoring iterations: 2Creating a new model with the removal of was_ket_ord and rec_pain_med due to insignificance.
glm_mod3 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Prop_dec_after_fent+Mode_CPOT+Con_Prop_Ben+Con_Prop_NMBA ,data=df_eda,family="gaussian")
summary(glm_mod3)##
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Prop_dec_after_fent +
## Mode_CPOT + Con_Prop_Ben + Con_Prop_NMBA, family = "gaussian",
## data = df_eda)
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.200 14.034 -0.371 0.71290
## Prop_dose_when_Dex10 38.119 30.742 1.240 0.22204
## Prop_dose_when_Dex100 31.309 22.416 1.397 0.17001
## Prop_dose_when_Dex21.2404 37.383 23.138 1.616 0.11384
## Prop_dose_when_Dex40 11.364 22.416 0.507 0.61491
## Prop_dose_when_Dex50 9.381 22.387 0.419 0.67739
## Prop_dose_when_Dex50.0002 16.529 22.387 0.738 0.46453
## Prop_dose_when_Dex51.2821 62.910 27.929 2.253 0.02971 *
## Prop_dose_when_Dex52 48.761 28.948 1.684 0.09969 .
## Prop_dose_when_Dex53.8642 60.659 22.387 2.710 0.00979 **
## Prop_dose_when_Dex54.1463 69.298 22.416 3.091 0.00357 **
## Prop_dose_when_Dex59.991 35.981 21.834 1.648 0.10700
## Prop_dose_when_Dex60 39.047 23.901 1.634 0.10998
## Prop_dose_when_Dex65 11.919 22.416 0.532 0.59778
## Prop_dose_when_Dex70 22.763 18.144 1.255 0.21675
## Prop_dose_when_Dex80 43.158 26.002 1.660 0.10459
## Prop_dose_when_Dex89.8261 29.158 21.834 1.335 0.18908
## Prop_dose_when_DexNotApp 20.831 14.252 1.462 0.15146
## Prop_dec_after_fent1 -20.944 15.054 -1.391 0.17165
## Prop_dec_after_fentNotApp -3.039 10.281 -0.296 0.76906
## Mode_CPOT1 -1.902 12.743 -0.149 0.88206
## Mode_CPOT2 1.866 11.114 0.168 0.86747
## Mode_CPOT3 9.024 9.122 0.989 0.32830
## Mode_CPOT4 10.521 10.082 1.044 0.30281
## Mode_CPOT5 29.029 8.541 3.399 0.00152 **
## Mode_CPOT6 31.886 12.971 2.458 0.01828 *
## Mode_CPOT7 NA NA NA NA
## Mode_CPOT8 2.301 18.826 0.122 0.90331
## Mode_CPOTNotApp -11.732 17.354 -0.676 0.50279
## Con_Prop_Ben -11.000 7.698 -1.429 0.16060
## Con_Prop_NMBA1 8.847 5.268 1.679 0.10068
## Con_Prop_NMBANotApp 6.785 22.387 0.303 0.76336
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 288.4435)
##
## Null deviance: 31735 on 71 degrees of freedom
## Residual deviance: 11826 on 41 degrees of freedom
## AIC: 635.63
##
## Number of Fisher Scoring iterations: 2Creating a new model with the removal of con_prop_ben and con_prop_NMBA due to insignificance.
glm_mod4 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Prop_dec_after_fent+Mode_CPOT ,data=df_eda,family="gaussian")
summary(glm_mod4)##
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Prop_dec_after_fent +
## Mode_CPOT, family = "gaussian", data = df_eda)
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6689 12.0325 0.056 0.95592
## Prop_dose_when_Dex10 35.6161 30.6364 1.163 0.25128
## Prop_dose_when_Dex100 27.5037 22.1297 1.243 0.22051
## Prop_dose_when_Dex21.2404 31.4240 22.1297 1.420 0.16266
## Prop_dose_when_Dex40 7.5586 22.1297 0.342 0.73431
## Prop_dose_when_Dex50 3.4673 20.4161 0.170 0.86592
## Prop_dose_when_Dex50.0002 21.5706 22.1297 0.975 0.33502
## Prop_dose_when_Dex51.2821 60.8641 27.8911 2.182 0.03448 *
## Prop_dose_when_Dex52 57.7160 27.8911 2.069 0.04442 *
## Prop_dose_when_Dex53.8642 65.7009 22.1297 2.969 0.00482 **
## Prop_dose_when_Dex54.1463 65.4926 22.1297 2.959 0.00495 **
## Prop_dose_when_Dex59.991 32.1317 20.9987 1.530 0.13313
## Prop_dose_when_Dex60 36.0523 23.7780 1.516 0.13662
## Prop_dose_when_Dex65 8.1141 22.1297 0.367 0.71563
## Prop_dose_when_Dex70 22.7545 17.5799 1.294 0.20230
## Prop_dose_when_Dex80 32.2657 25.1141 1.285 0.20560
## Prop_dose_when_Dex89.8261 25.3087 20.9987 1.205 0.23455
## Prop_dose_when_DexNotApp 18.7854 13.7528 1.366 0.17890
## Prop_dec_after_fent1 -22.5511 14.5088 -1.554 0.12727
## Prop_dec_after_fentNotApp -5.1024 9.9406 -0.513 0.61031
## Mode_CPOT1 -4.7385 12.5270 -0.378 0.70706
## Mode_CPOT2 -1.3287 10.3857 -0.128 0.89878
## Mode_CPOT3 8.9382 9.1286 0.979 0.33286
## Mode_CPOT4 10.5658 9.6682 1.093 0.28041
## Mode_CPOT5 28.2192 8.5576 3.298 0.00194 **
## Mode_CPOT6 38.9734 12.5270 3.111 0.00327 **
## Mode_CPOT7 NA NA NA NA
## Mode_CPOT8 -10.4584 17.4390 -0.600 0.55178
## Mode_CPOTNotApp -13.4918 17.4390 -0.774 0.44327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 294.3858)
##
## Null deviance: 31735 on 71 degrees of freedom
## Residual deviance: 12953 on 44 degrees of freedom
## AIC: 636.18
##
## Number of Fisher Scoring iterations: 2Performing variable selection using deviance tables
anova(glm_mod4,test="Chisq")## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Per_prop_over50
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 71 31735
## Prop_dose_when_Dex 17 10969.2 54 20766 0.003101 **
## Prop_dec_after_fent 2 1061.9 52 19704 0.164695
## Mode_CPOT 8 6750.8 44 12953 0.003453 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The final model for the % of time of propofol above 50 mcg/kg/min (Per_prop_over50) based on the deviance table can be written as follows:
Per_prop_over50 = 0.241 + 2.961 * Prop_dose_when_Dex + 3.334 * Mode_CPOT
glm_mod5 <- glm(Per_prop_over50~ Prop_dose_when_Dex+Mode_CPOT,data=df_eda,family="gaussian")
summary(glm_mod5)##
## Call:
## glm(formula = Per_prop_over50 ~ Prop_dose_when_Dex + Mode_CPOT,
## family = "gaussian", data = df_eda)
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9746 11.8023 0.167 0.86786
## Prop_dose_when_Dex10 10.8880 27.1926 0.400 0.69071
## Prop_dose_when_Dex100 21.0956 20.9709 1.006 0.31970
## Prop_dose_when_Dex21.2404 25.0159 20.9709 1.193 0.23903
## Prop_dose_when_Dex40 1.1504 20.9709 0.055 0.95649
## Prop_dose_when_Dex50 6.6714 20.2591 0.329 0.74342
## Prop_dose_when_Dex50.0002 15.1625 20.9709 0.723 0.47332
## Prop_dose_when_Dex51.2821 53.5848 27.1926 1.971 0.05481 .
## Prop_dose_when_Dex52 50.4366 27.1926 1.855 0.07004 .
## Prop_dose_when_Dex53.8642 59.2928 20.9709 2.827 0.00693 **
## Prop_dose_when_Dex54.1463 59.0844 20.9709 2.817 0.00711 **
## Prop_dose_when_Dex59.991 30.2333 20.2591 1.492 0.14244
## Prop_dose_when_Dex60 28.3028 22.4704 1.260 0.21418
## Prop_dose_when_Dex65 1.7060 20.9709 0.081 0.93552
## Prop_dose_when_Dex70 18.2269 17.4245 1.046 0.30100
## Prop_dose_when_Dex80 24.9863 24.2733 1.029 0.30869
## Prop_dose_when_Dex89.8261 23.4104 20.2591 1.156 0.25383
## Prop_dose_when_DexNotApp 11.5061 11.7672 0.978 0.33328
## Mode_CPOT1 -3.8672 12.6234 -0.306 0.76072
## Mode_CPOT2 -0.4575 10.4532 -0.044 0.96528
## Mode_CPOT3 9.8094 9.1777 1.069 0.29072
## Mode_CPOT4 6.0560 9.3834 0.645 0.52188
## Mode_CPOT5 29.5605 8.6006 3.437 0.00126 **
## Mode_CPOT6 39.8446 12.6234 3.156 0.00282 **
## Mode_CPOT7 NA NA NA NA
## Mode_CPOT8 -9.5872 17.5953 -0.545 0.58847
## Mode_CPOTNotApp -12.6206 17.5953 -0.717 0.47683
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 300.4846)
##
## Null deviance: 31735 on 71 degrees of freedom
## Residual deviance: 13822 on 46 degrees of freedom
## AIC: 636.86
##
## Number of Fisher Scoring iterations: 2Running diagnostics and checking model adequacy
par(mfrow=c(2,2))
plot(glm_mod5)## Warning: not plotting observations with leverage one:
## 1, 2, 4, 7, 8, 9, 12, 17, 18, 19, 20, 21, 22, 26, 71## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
Evaluating evidence of interaction between a quantitative and a qualitative variable
There is evidence of interaction between a quantitative predictor and a qualitative predictor based off the EDA. There seems to be some evidence of interaction between the total days on propofol infusion (Total_days_prop) (quantitative) and the organization (Org_Name) (qualitative). An ANCOVA will be run.
#Starting with ANOVA
Propdata_anova = select(df_eda, Per_prop_over50,Org_Name, Total_days_prop)
glm_anova <- glm(Per_prop_over50~ factor(Org_Name),data=Propdata_anova,family="gaussian")
summary(glm_anova)##
## Call:
## glm(formula = Per_prop_over50 ~ factor(Org_Name), family = "gaussian",
## data = Propdata_anova)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.047 9.426 1.172 0.245
## factor(Org_Name)IU HBloomH 3.073 11.329 0.271 0.787
## factor(Org_Name)IU HBMH 19.147 9.886 1.937 0.057 .
## factor(Org_Name)IU HMH -4.314 10.539 -0.409 0.684
## factor(Org_Name)IU HUH -2.468 14.399 -0.171 0.864
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 355.4004)
##
## Null deviance: 31735 on 71 degrees of freedom
## Residual deviance: 23812 on 67 degrees of freedom
## AIC: 634.02
##
## Number of Fisher Scoring iterations: 2There is a small amount of significance when the organization name is ‘IU HBMH’, which matches with the findings so far.
#Running ANCOVA
glm_ancova <- glm(Per_prop_over50~ factor(Org_Name) + Total_days_prop, data=Propdata_anova,family="gaussian")
summary(glm_ancova)##
## Call:
## glm(formula = Per_prop_over50 ~ factor(Org_Name) + Total_days_prop,
## family = "gaussian", data = Propdata_anova)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.64986 10.11925 1.151 0.2538
## factor(Org_Name)IU HBloomH 2.98745 11.42234 0.262 0.7945
## factor(Org_Name)IU HBMH 19.00153 9.99409 1.901 0.0616 .
## factor(Org_Name)IU HMH -4.36904 10.62054 -0.411 0.6821
## factor(Org_Name)IU HUH -2.65121 14.54262 -0.182 0.8559
## Total_days_prop -0.05482 0.31810 -0.172 0.8637
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 360.623)
##
## Null deviance: 31735 on 71 degrees of freedom
## Residual deviance: 23801 on 66 degrees of freedom
## AIC: 635.99
##
## Number of Fisher Scoring iterations: 2anova(glm_ancova)## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: Per_prop_over50
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 71 31735
## factor(Org_Name) 4 7923.1 67 23812
## Total_days_prop 1 10.7 66 23801Once again, there is a small amount of significance when the organization name is ‘IU HBMH’, which matches with the findings so far.
#interaction between the total days on propofol infusion (Total_days_prop) (quantitative) and the organization (Org_Name) (qualitative)
par(mar=c(1,1,1,1))
par(mfrow=c(1,1))
ggplot(data = Propdata_anova, aes(x=Total_days_prop, y=Per_prop_over50, color=factor(Org_Name)))+
geom_point() + geom_smooth() + xlab("total days on propofol infusion (Total_days_prop)") + ylab("% of time of propofol above 50 mcg/kg/min (Per_prop_over50)")## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 4.93## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 8.07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 145.68## Warning in sqrt(sum.squares/one.delta): NaNs produced## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 4.93## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 8.07## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 145.68## Warning in stats::qt(level/2 + 0.5, pred$df): NaNs produced## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 4## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1.1299e-17## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 121## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at 4## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 1.1299e-17## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 121## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 3.97## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 5.03## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 1.0609## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : span too small. fewer
## data values than degrees of freedom.## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 3.97## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 5.03## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 0## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 1.0609## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
#testing for interaction
glm_ancova1 <- glm(Per_prop_over50~ factor(Org_Name) + Total_days_prop +factor(Org_Name)*Total_days_prop, data=Propdata_anova, family="gaussian")
summary(glm_ancova1)##
## Call:
## glm(formula = Per_prop_over50 ~ factor(Org_Name) + Total_days_prop +
## factor(Org_Name) * Total_days_prop, family = "gaussian",
## data = Propdata_anova)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.9435 21.5036 0.974 0.334
## factor(Org_Name)IU HBloomH 2.7782 23.5070 0.118 0.906
## factor(Org_Name)IU HBMH 6.6013 21.9946 0.300 0.765
## factor(Org_Name)IU HMH -12.6945 23.2714 -0.545 0.587
## factor(Org_Name)IU HUH -19.5092 40.3735 -0.483 0.631
## Total_days_prop -0.8997 1.7499 -0.514 0.609
## factor(Org_Name)IU HBloomH:Total_days_prop -0.1170 1.9014 -0.062 0.951
## factor(Org_Name)IU HBMH:Total_days_prop 1.2169 1.7991 0.676 0.501
## factor(Org_Name)IU HMH:Total_days_prop 0.7481 1.9037 0.393 0.696
## factor(Org_Name)IU HUH:Total_days_prop 1.8315 4.5654 0.401 0.690
##
## (Dispersion parameter for gaussian family taken to be 367.4723)
##
## Null deviance: 31735 on 71 degrees of freedom
## Residual deviance: 22783 on 62 degrees of freedom
## AIC: 640.84
##
## Number of Fisher Scoring iterations: 2#running model diagnostics
anova(glm_ancova, glm_ancova1)## Analysis of Deviance Table
##
## Model 1: Per_prop_over50 ~ factor(Org_Name) + Total_days_prop
## Model 2: Per_prop_over50 ~ factor(Org_Name) + Total_days_prop + factor(Org_Name) *
## Total_days_prop
## Resid. Df Resid. Dev Df Deviance
## 1 66 23801
## 2 62 22783 4 1017.8The interaction term is significant as the difference in deviance between the larger model and the smaller model is statistically significant.
Final model with interaction term is as follows:
Per_prop_over50 = 0.241 + 2.961 * Prop_dose_when_Dex + 3.334 * Mode_CPOT - 0.518 * Total_days_prop
The final explanatory variables for the % of time of propofol above 50 mcg/kg/min are Prop_dose_when_dex, Mode_CPOT, and Total_days_prop.
Bibliography
Physiopedia. (n.d.). Richmond agitation-sedation scale (RASS). Physiopedia. https://www.physio-pedia.com/Richmond_Agitation-Sedation_Scale_(RASS) Description of Richmond Agitation-Sedation Scale (RASS).
South Tees Hospitals NHS FT. (2023, January 20). CPOT = Critical-Care Pain Observation Tool. Critical Care Services - The James Cook University Hospital. https://ccs-sth.org/resources/Documents/Sedation Analgesia Delirium in CC/CCS STH Guideline template CPOT.pdf Description of CPOT = Critical-Care Pain Observation Tool.